Práctica II

Análisis clúster

Instrucciones (leer antes de empezar)

• Modifica dentro del documento .Rmd tus datos personales (nombre y DNI) ubicados en la cabecera del archivo.

• Asegúrate antes de seguir editando el documento que el archivo .Rmd compila (Knit) correctamente y se genera el html correspondiente.

• Los chunks creados están o vacíos o incompletos, de ahí que tengan la opción eval = FALSE. Una vez que edites lo que consideres debes de cambiar a eval = TRUE para que los chunk se ejecuten

Paquetes necesarios

Necesitaremos los siguientes paquetes:

• Manejo de datos: paquete {tidyverse}.
• Modelos: paquete {tidymodels}
• Lectura excel: paquete {readxl}
• Resumen numérico: paquete {skimr}.
• Visualización de clústers y PCA: paquete {factoextra} y {FactoMineR}
• Clustering divisivo: paquete {cluster}

# Borramos variables del environment
rm(list = ls())
library(tidyverse)
library(tidymodels)
library(readxl)
library(skimr)
library(factoextra)
library(FactoMineR)
library(cluster)

library(corrplot)
library(scales)
library(heatmaply)

Carga de datos

El archivo de datos a usar será provincias.xlsx

provincias <- read_xlsx(path = "./provincias.xlsx")

El fichero contiene información socioeconómica de las provincias españolas

glimpse(provincias)

Rows: 52
Columns: 19
$ Prov          <chr> "Albacete", "Alicante", "Almería", "Álava", "A…
$ Poblacion     <dbl> 396987, 1868438, 701688, 321932, 1061756, 6909…
$ Mortalidad    <dbl> 9.41, 8.00, 6.91, 7.72, 12.16, 9.54, 7.08, 8.1…
$ Natalidad     <dbl> 9.02, 8.52, 11.41, 10.26, 6.26, 8.96, 9.51, 9.…
$ IPC           <dbl> 101.819, 102.309, 101.454, 102.743, 102.070, 1…
$ NumEmpresas   <dbl> 26701, 130438, 40327, 19548, 67451, 39640, 871…
$ Industria     <dbl> 3332, 9992, 2191, 2048, 3496, 3059, 4413, 2741…
$ Construccion  <dbl> 3428, 17418, 5388, 2669, 8435, 4599, 14485, 50…
$ CTH           <dbl> 11111, 52720, 18184, 7557, 28045, 18420, 30644…
$ Infor         <dbl> 189, 1851, 366, 323, 787, 332, 1403, 11134, 15…
$ AFS           <dbl> 624, 2745, 944, 354, 1445, 948, 1627, 8829, 16…
$ APT           <dbl> 3278, 19865, 5840, 3341, 10598, 5235, 16444, 8…
$ TasaActividad <dbl> 57.72, 58.21, 60.63, 57.86, 51.42, 56.00, 68.6…
$ TasaParo      <dbl> 23.28, 23.39, 31.08, 13.33, 16.97, 29.63, 13.8…
$ Ocupados      <dbl> 144.8, 687.1, 234.4, 133.0, 389.0, 224.6, 555.…
$ PIB           <dbl> 7235991, 32139347, 11909684, 10716021, 2177043…
$ CANE          <dbl> 20888, 25968, 22945, 3691, 23910, 37438, 10748…
$ TVF           <dbl> 215948, 1274096, 395086, 155767, 613905, 37249…
$ VS            <dbl> 30244, 326705, 72486, 9791, 73250, 49441, 8571…

Algunas de las variables son:

• Prov: nombre de la provincia
• Poblacion: habitantes
• Mortalidad, Natalidad: tasa de mortalidad/natalidad (en tantos por mil)
• IPC: índice de precios de consumo (sobre un valor base de 100).
• NumEmpresas: número de empresas.
• PIB: producto interior bruto
• CTH: coyuntura turística hotelera (pernoctaciones en establecimientos hoteleros)

Ejercicio 1:

Calcula la matriz de covarianzas y de correlaciones. Usa el paquete {corrplot} para una representación gráfica de la misma. Detalla y comenta lo que consideres para su correcta interpretación.

cov_mat <- provincias %>% select(-Prov) %>% cov() # quitamos la variab (id) categórica

cor_mat <- provincias %>% select(-Prov) %>% cor()

head(cor_mat[c(1,2,3), c(1,2,3)])

            Poblacion Mortalidad  Natalidad
Poblacion   1.0000000 -0.3435994  0.1144594
Mortalidad -0.3435994  1.0000000 -0.7365292
Natalidad   0.1144594 -0.7365292  1.0000000

corrplot(cor_mat, type = "upper",
         tl.col = "black",  method = "ellipse")

# Para interpretar este gráfico, recordemos que se trata de la representación de la
# matriz de correlaciones entre las variables numéricas de nuestro dataset. La
# matriz es, lógicamente, cuadrada (18x18 en nuestro caso) y simétrica, y realmente
# nos vasta con representar solo uno de los lados de la diagonal principal para tener
# todas las correlaciones representadas. Obviamente, una variable consigo misma 
# tiene una corr = 1 exacta, por lo que los valores en la diagonal corresponden
# a corr = 1, y por eso se pintan con el tono de azul correspondiente a 1, y como 
# una línea recta. El resto de corr se representan también utilizando tanto el 
# color como la forma de elipse, con sus ejes más desiguales cuanto mayor es la
# corr (hasta que corr = 1 corresponde a eje menor = 0, ie una línea recta).
# Además, para indicar si la corr es positiva o negativa, además del color, se juega 
# con la orientación de la elipse.
# 
# En nuestro gráfico, vemos que hay varias variables altamente correlacionadas entre sí,
# todas ellas positivamente, y que, salvo por TVF, todas las variabs altamente corr son de 
# carácter económico (el dataset tiene variables tanto económicas como poblacionales).
# 
# Por otra parte, se ven también bastantes correlaciones bajas, tanto positivas como
# negativas. Y corr altas negativas hay muy pocas: entre Natalidad y Mortalidad, entre
# TasaActividad y Mortalidad, y entre IPC y TasaParo.

Ejercicio 2:

Calcula con eigen() los autovalores y autovectores de la matriz de correlaciones e interpreta dichos resultados en relación a las componentes principales de las variables originales.

autocosas <- eigen(cor_mat) # autovalores y autovectores

autocosas$values # autovalores

 [1] 1.146640e+01 2.560508e+00 1.634097e+00 9.340488e-01 4.565320e-01
 [6] 4.143630e-01 3.071743e-01 1.166154e-01 7.309913e-02 2.034496e-02
[11] 9.027568e-03 3.566115e-03 2.372720e-03 8.013952e-04 5.472749e-04
[16] 3.512710e-04 1.028975e-04 4.950858e-05

autocosas$vectors # autovectores

             [,1]         [,2]         [,3]        [,4]         [,5]
 [1,]  0.29352698  0.002493342  0.049947340 -0.05339540 -0.015898841
 [2,] -0.10629199 -0.527120141  0.188867032 -0.16131678  0.024146585
 [3,]  0.04062702  0.495419403 -0.270681171 -0.10967036 -0.382061614
 [4,]  0.10991169 -0.365322291 -0.262317619  0.43504738 -0.556757431
 [5,]  0.29418239 -0.026147510  0.007952241 -0.06913798 -0.006367061
 [6,]  0.28562283 -0.044737470  0.046369149  0.02345108 -0.150042973
 [7,]  0.29326489 -0.045381770 -0.012082096 -0.02639413  0.008113738
 [8,]  0.29286608 -0.010617503  0.048841149 -0.02758090 -0.028366126
 [9,]  0.28150192 -0.042018125 -0.064640257 -0.22152688  0.070433510
[10,]  0.29246592 -0.016454205  0.039554303 -0.09168100  0.003041574
[11,]  0.29072490 -0.029433207 -0.028068230 -0.14186101  0.031732816
[12,]  0.11433729  0.330614591 -0.362734482  0.46305159  0.133043722
[13,] -0.01399634  0.461562171  0.387356995 -0.22020186 -0.055362718
[14,]  0.29442838 -0.016838806  0.002429338 -0.06005252 -0.001232463
[15,]  0.29090708 -0.036157691 -0.037491740 -0.13360720 -0.008622583
[16,]  0.01778860  0.096184477  0.656743135  0.27839359 -0.488037746
[17,]  0.29156666 -0.001599604  0.099761847  0.04414466  0.046676773
[18,]  0.17234327  0.047821027  0.290142940  0.56712077  0.502605369
              [,6]         [,7]         [,8]         [,9]
 [1,]  0.009261485 -0.067375388  0.049695939  0.063533937
 [2,] -0.033735542  0.194593899 -0.677240259  0.377641454
 [3,]  0.313460111  0.588024487 -0.256257868  0.063132926
 [4,]  0.350737255 -0.328406660 -0.146407863 -0.180898036
 [5,] -0.006803308 -0.016541352  0.004345160  0.087585774
 [6,]  0.061617223  0.064730810  0.277225726  0.598423544
 [7,] -0.006262106  0.026640155  0.030059205  0.029546241
 [8,]  0.047239229 -0.068140824  0.107571562  0.261268280
 [9,] -0.125617691  0.051367829 -0.249587279 -0.415131011
[10,] -0.054714373  0.018516638 -0.071851535 -0.178586941
[11,] -0.080055846  0.008053971 -0.124821218 -0.199015131
[12,] -0.528510461 -0.186075273 -0.370767601  0.246498807
[13,]  0.304199666 -0.591402972 -0.343449492  0.089255270
[14,] -0.039890164 -0.023453829  0.005640881 -0.009548811
[15,] -0.077030649  0.008952134 -0.056398210 -0.156415106
[16,] -0.419227586  0.199631828  0.013485836 -0.149563592
[17,]  0.073722587 -0.005584984  0.069914036  0.095085695
[18,]  0.430852254  0.257315614 -0.125215945 -0.114868372
             [,10]        [,11]         [,12]        [,13]
 [1,]  0.350085935  0.186230539  0.1711785775  0.010288303
 [2,]  0.075589860 -0.007930229  0.0420207640  0.021526962
 [3,]  0.082446318 -0.020289162  0.0162179700 -0.017268779
 [4,] -0.024453594  0.018055215 -0.0165433445 -0.009256375
 [5,]  0.022328838 -0.155683246 -0.1868313262 -0.295730236
 [6,] -0.549771328  0.313257119 -0.0983659710  0.159247816
 [7,] -0.191834125 -0.821885503  0.1664388522  0.342872601
 [8,]  0.275273747 -0.176716870 -0.2371619940 -0.377496712
 [9,] -0.455403694  0.145292814 -0.1647985653 -0.169331304
[10,]  0.280705914  0.197698905 -0.4202117212  0.708492865
[11,] -0.063053183 -0.018374368 -0.2760928653 -0.285162833
[12,]  0.002558714  0.022324556  0.0006653841  0.022821245
[13,] -0.112317487 -0.013947370  0.0209740158  0.049593889
[14,]  0.154687198 -0.048862657  0.2063706502 -0.018496220
[15,] -0.098205180  0.224212469  0.7033944369  0.025058312
[16,] -0.018061942 -0.036012420  0.0068803352 -0.057390515
[17,]  0.320576240  0.138218764  0.1485323831 -0.081502147
[18,] -0.097744721  0.033767221  0.0002578474 -0.018041259
             [,14]        [,15]         [,16]         [,17]
 [1,]  0.174252449  0.528193532 -0.5074390165 -0.3645836076
 [2,]  0.012714926  0.006179773 -0.0007385454 -0.0017830808
 [3,]  0.007145884  0.010947891 -0.0010884443  0.0031311475
 [4,]  0.003959246  0.003740459  0.0067705625  0.0041585701
 [5,] -0.198427286 -0.216945648  0.1011828614 -0.2158364244
 [6,]  0.083661952 -0.036645845 -0.0274895935 -0.0289281260
 [7,]  0.040333104  0.172074291  0.0517103111 -0.1011674569
 [8,] -0.459614346  0.104642291 -0.1348807749  0.3944508893
 [9,]  0.063539475  0.453142059  0.0509083701  0.3227546858
[10,] -0.207133987 -0.137066568  0.0335887182  0.0637126478
[11,]  0.214532746 -0.367229939 -0.0374738932 -0.5220390566
[12,] -0.007663637  0.017628269  0.0166789316 -0.0042896852
[13,]  0.006561163 -0.031150623  0.0136695061 -0.0010611951
[14,]  0.550167465 -0.425843706 -0.2681093335  0.5192283126
[15,] -0.491058860 -0.230864767  0.0185550861 -0.0582894917
[16,] -0.009889913  0.008894585 -0.0107381109  0.0036300393
[17,]  0.265988870  0.187180605  0.7868348317 -0.0040644103
[18,] -0.049235072 -0.043019092 -0.1164353064  0.0006033047
              [,18]
 [1,] -0.1033819874
 [2,]  0.0093642618
 [3,]  0.0024522685
 [4,]  0.0006333815
 [5,] -0.7807168763
 [6,]  0.0538942552
 [7,]  0.0824471845
 [8,]  0.3488906784
 [9,] -0.0852073738
[10,] -0.0228329277
[11,]  0.4629930678
[12,]  0.0075771232
[13,]  0.0059635657
[14,] -0.1113370996
[15,]  0.0871095987
[16,] -0.0033788830
[17,]  0.0775931859
[18,] -0.0043631422

# El dataset original contiene 19 variabs, pero una es la provincia (no numérica),
# por lo que la matriz de correlaciones es 18x18.

# - Los autovectores (loadings) son las coordenadas de las componentes principales
#   respecto a las variabs orig.
#   Cuanto mayor es el autovalor asociado a cada autovector, mayor es la información
#   (varianza en los datos) que captura ese autovector (esa dirección en el espacio
#   de variables numéricas estudiado).

# - No hay autovalores (ni por ende, autovectores) repetidos, lo que significa que
#   todas las variabs orig contienen cierta información que no está capturada en las
#   demás variabs. Así, podemos tener hasta 18 nuevas comps principales.
#   Podemos obtener la "varianza explicada acumulada" por cada componente así:

cumsum_var <- cumsum(autocosas$values) / sum(autocosas$values)
cumsum_var # 6 comps explican el 97% de la varianza

 [1] 0.6370222 0.7792726 0.8700558 0.9219474 0.9473102 0.9703304
 [7] 0.9873957 0.9938743 0.9979353 0.9990656 0.9995672 0.9997653
[13] 0.9998971 0.9999416 0.9999720 0.9999915 0.9999972 1.0000000

#   Con esto, vemos que con solo 6 comps principales ya explicamos más del 95%
#   de la varianza.

Ejercicio 3:

Haciendo uso de PCA() del paquete {FactoMineR} calcula todas las componentes principales. Repite de nuevo el análisis con el mínimo número de componentes necesarias para capturar al menos el 95% de la información de los datos.

# calculamos todas las comp princ:
pca_fit <-  PCA(provincias[, -1],
                scale.unit = TRUE, # pues los datos orig no están estandarizados
                ncp = 18,
                graph = FALSE)

# buscamos las que explican el 95% de la varianza:
pca_fit$eig[ ,3] # varianza acumulada --> con 5 comp llegamos al 94.7%; con 6, al 97.03%

   comp 1    comp 2    comp 3    comp 4    comp 5    comp 6    comp 7 
 63.70222  77.92726  87.00558  92.19474  94.73102  97.03304  98.73957 
   comp 8    comp 9   comp 10   comp 11   comp 12   comp 13   comp 14 
 99.38743  99.79353  99.90656  99.95672  99.97653  99.98971  99.99416 
  comp 15   comp 16   comp 17   comp 18 
 99.99720  99.99915  99.99972 100.00000 

# repetimos el PCA, pero con ncp = 6
pca_fit <-  PCA(provincias[, -1],
                scale.unit = TRUE,
                ncp = 6,
                graph = FALSE)

Ejercicio 4:

Realiza las gráficas y análisis numéricos que consideres más útiles para poder interpretar adecuadamente al análisis de componentes principales obtenidas. Detalla todo lo que puedas concluir y, en concreto, responde a las siguientes preguntas

# 1) varianza explicada por cada comp:
fviz_eig(pca_fit,
         choice = "variance",
         ncp = 6, # nº de autovalores a pintar
         barfill = "darkolivegreen",
         addlabels = TRUE) +
  theme_minimal() +
  labs() +
  xlab("Componentes") +
  ylab("% de varianza explicada") +
  ggtitle("Varianza explicada por componentes")

# Vemos que con tan sólo una comp principal ya capturamos casi el 64% de la varianza
# de los datos originales.
# El 36% restante, no obstante, está más repartido entre el resto de comps.

# 2) varianza explicada ACUMULADA:
cumvar <- as_tibble(pca_fit$eig)
names(cumvar) <- c("lambda", "var", "cumvar")

ggplot(cumvar[1:6, ],
       aes(x=1:6, y=cumvar)) +
  geom_col(fill = "pink") +
  geom_hline(yintercept = 95) + # pintar el 95%
  theme_minimal() +
  scale_x_discrete(limit = c("Dim 1","Dim 2","Dim 3","Dim 4","Dim 5","Dim 6")) +
  labs() +
  xlab("Componentes") +
  ylab("% de varianza ACUMULADA explicada") +
  ggtitle("Varianza ACUMULADA explicada por componentes")

# Este gráfico contiene información muy similar al anterior, con la diferencia de
# que nos permite visualizar directamente cuántas comps principales hacen falta para
# superar el umbral del 95% de varianza acumulada explicada.

# 3.a) proyección bidim de las variabs orig sobre las comp 1 y 2 (las de mayor var explicada)
col <- c("#00AFBB", "#E7B800", "#FC4E07")

fviz_pca_var(pca_fit, col.var = "cos2", # por dfto axes = 1:2 (comps 1 y 2)
             gradient.cols = col,
             repel = TRUE) +
  theme_minimal() + 
  labs(color= "Peso") +
  ggtitle("Coordenadas de las variables en las dos Comp Principales")

# Lo primero que destaca es el hecho de que, de las 18 variabs orig, 11 están
# altamente correlacionadas con la comp 1. Esto concuerda con lo que vimos
# anteriormente de que solamente la primera comp ya capturaba casi el 64% de
# la varianza. También tiene sentido que esto sea así, ya que, si recordamos
# el corrplot de la matriz de correlaciones cor_mat, veíamos que muchas de
# las variabs orig estaban altamente correlacionadas entre sí, y por tanto
# tiene sentido que una sola comp principal pueda explicarlas a todas.

# Por su parte, las otras 7 variabs muestran comportamientos diferentes, y parecen
# estar relativamente bien explicadas por la combinación de Dim.1 y Dim.2 (recordemos
# que estas dos comps juntas suman más del 75% de la varianza acumulada explicada), con
# la clara salvedad de la variab CANE.

# Por último, recordemos que nuestro dataset tiene variabs económicas y poblacionales.
# Entonces, cabe destacar que las 11 variabs correladas con la comp 1 son todas de
# carácter económico, a excepción de "Poblacion".
# Las demás variabs económicas, que no están altamente correlacionadas con la comp 1,
# son: IPC, VS, TasaParo y TasaActividad.
# Esta información nos será útil más adelante.

# 3.b) 
# Podemos buscar en qué comps están mejor representadas estas variabs. Por ej, CANE:

which(colnames(provincias[,-1]) == "CANE") # identificar CANE (sin tener en cuenta

[1] 16

# el número col de la variab categórica, pues no aparece luego en pca_fit)

pca_fit$var$cos2[16, ] # buscamos en qué comps está mejor representada la variab CANE

      Dim.1       Dim.2       Dim.3       Dim.4       Dim.5 
0.003628361 0.023688421 0.704804777 0.072391578 0.108737179 
      Dim.6 
0.072825031 

as.data.frame(pca_fit$var$cos2[16, ]) %>%
  slice_max(pca_fit$var$cos2[16, ],
            n = 2) # resultado --> Dim.3 y Dim.5

      pca_fit$var$cos2[16, ]
Dim.3              0.7048048
Dim.5              0.1087372

col <- c("#00AFBB", "#E7B800", "#FC4E07")
fviz_pca_var(pca_fit, col.var = "cos2",
             axes = c(3,5),
             gradient.cols = col,
             repel = TRUE) +
  theme_minimal() + 
  labs(color= "Peso") +
  ggtitle("Coordenadas de las variables en Dim.3 y Dim.5")

# Vemos cómo Dim.3 y Dim.5 capturan mejor la info de la variab CANE

### (De forma análoga, podríamos estudiar el resto de variabs y comps)

¿Cuál es la expresión para calcular la primera componente en función de las variables originales?

# los autovectores (loadings) son las coordenadas de las componentes respecto a las variabs orig

pca_fit$svd$V[ ,1] # coords de la comp 1

 [1]  0.29352698 -0.10629199  0.04062702  0.10991169  0.29418239
 [6]  0.28562283  0.29326489  0.29286608  0.28150192  0.29246592
[11]  0.29072490  0.11433729 -0.01399634  0.29442838  0.29090708
[16]  0.01778860  0.29156666  0.17234327

\[Dim_1 = 0.296*Poblacion - 0.106*Mortalidad + 0.041*Natalidad + 0.110*IPC + 0.294*NumEmpresas + 0.286*Industria + 0.293*Construccion + 0.293*CTH + 0.282*Infor + 0.292*AFS + 0.291*APT + 0.114*TasaActividad - 0.014*TasaParo + 0.294*Ocupados + 0.291*PIB + 0.018*CANE + 0.292*TVF + 0.172*VS \]

Proporciona las nuevas coordenadas de los datos.

# las nuevas coordenadas (o scores) de los datos en la base de 6 comps principales es:
pca_fit$ind$coord

         Dim.1         Dim.2       Dim.3       Dim.4       Dim.5
1  -1.41017686  0.4734423651  0.09617150 -0.45796713  0.13095494
2   3.38367603  0.5397505581  1.91865320  2.42036146 -1.93760003
3  -0.61651226  2.6139666809  0.20772292 -0.14164444 -0.04194355
4  -1.44429594 -0.0008444871 -2.03183179 -0.11310052  0.46785078
5  -0.20350137 -1.9525986131  1.29786007 -0.71388934 -0.36860431
6  -1.04846596  1.1677693270  1.79554023 -0.85440440  0.07544289
7   1.52621574  0.2601071540 -2.51922849  2.36354690  0.48085217
8  13.68268794 -1.6121155618 -0.86679541 -0.27855369  1.03981506
9   0.57644373 -1.5081871153 -1.17979880 -0.38255980  0.79824422
10 -1.20182455 -1.5496195139 -0.89154966  0.98854013  0.60182044
11 -0.84908539 -1.1259429948 -0.59433092  0.40115869 -0.13074923
12 -0.68987702  0.8209167732  0.51814851  0.38168516 -0.54791911
13 -2.12522644  3.3258586141 -1.81108385 -1.16466618  0.29805470
14 -1.39218950  0.8151885534  1.81886692 -0.77389817  0.21646631
15  0.63521316 -1.4418586784  0.75867607  0.24388244  0.25788234
16 -2.13228382 -0.5693293626  1.04829137 -0.93167013 -0.47752664
17 -1.50273216 -0.1797543487  1.26924013 -0.22326789 -0.04115631
18  0.12808408  1.7712185679  0.36869418 -0.24856268 -0.25948344
19 -0.59357443  0.8106412847  1.29151505 -0.16892787  1.08430886
20 -0.28110176 -1.4851526175 -1.87915578  0.26633383  1.18449911
21  0.38750851  0.5444284560 -1.38741262  1.63782864 -0.62910775
22 -0.07240306  1.1235420649  1.51103758  0.50031682  0.72301216
23 -1.40758889  1.5419225448 -1.84940927  0.90565227 -0.78769934
24 -1.26474244  1.1675386017  0.01287591 -0.32743518  0.07250879
25 -1.77555839 -0.9839009562 -0.58738677  0.30393125 -0.29036951
26 -1.22127859  1.4243451642  3.40732381 -0.47898046  1.42770039
27 -1.46356388 -2.0160745559  0.87729273 -0.80858387 -0.23952614
28 -0.83527088 -0.1355632650 -1.47185005  1.18507362  0.74051659
29 -1.82595354 -2.7845330165  0.87126549 -0.35278099  0.32404978
30 16.77829847 -0.3659310780 -0.84914740 -2.36482726 -0.88049226
31 -2.21794591  4.7821956680 -1.90516137 -2.23082541  0.71383090
32  1.52219165  1.4419889136  0.58016589  0.90077234  0.43650450
33  2.00578012  1.3246290015  0.86903425  1.10396096 -0.35663749
34 -0.65330322  0.0783395528 -1.09566702 -0.08112406  0.22905253
35 -1.96532314 -2.8583465650  1.09835557 -1.20421552 -0.32532747
36 -2.12245758 -1.9505762018 -0.69512991 -0.14223589  0.04978965
37  0.09184283  1.8568628084 -0.33025189 -0.82364779 -1.37494830
38  0.03605515 -0.6068176749  0.05204844 -0.32771414 -0.01608646
39 -1.38317512 -0.4839449757 -1.35401168  0.28829039  0.20077731
40 -1.61199535 -1.4247515635 -0.12106972 -0.08568660 -0.57711355
41 -0.02919167  1.5729750162 -0.12649969 -0.43206728 -1.55965419
42 -1.93110326 -1.0149596929 -1.10961470  0.26282925 -0.26624987
43  1.94829323  1.7745711703  0.71196983 -0.54553119  0.54866804
44 -2.39905891 -1.8572255862 -0.77794225 -0.50345472 -0.65632722
45  0.17541169  1.0403567836 -0.10195495  1.51207804 -0.53700307
46 -2.18541694 -1.0816086840 -0.35121773 -0.41279186 -0.51591398
47 -0.46148674  1.4488771551  0.81164464  0.46323028  0.04743993
48  4.76983184  0.3601408245  2.96073023  2.34944835  0.91936057
49 -1.00710369 -0.7041391017 -1.05224187  0.19559460  0.29014683
50 -2.28666239 -3.1694811349  0.57775948 -0.62184293  0.16723585
51  0.11540844 -0.2367335810 -0.15614346  0.06774386  0.14172596
52 -2.15151156 -0.9815826773  0.36500304 -0.54540190 -0.85107240
         Dim.6
1  -0.37892322
2   1.79860288
3   0.05037528
4   0.09437501
5   0.01671182
6  -0.57465897
7  -0.56722730
8   0.76543430
9   0.46557476
10  0.53211012
11  0.39106047
12 -0.29386801
13  0.64775615
14 -0.56462157
15 -0.04619127
16 -0.08047104
17  0.07416747
18  1.27315904
19 -0.01365680
20  0.48379040
21  0.50800878
22 -0.08984337
23 -1.01677677
24  0.92998978
25 -0.79299337
26 -1.00010051
27  0.66516503
28 -0.92356915
29 -0.16130804
30 -0.97422963
31  1.40966039
32  0.40629457
33  0.73782389
34 -0.70374812
35  0.88838418
36  0.04140390
37 -0.54986198
38  0.41943263
39 -0.24302453
40 -0.24339861
41 -0.74670109
42 -0.39916228
43 -0.20509500
44 -0.41982505
45 -0.22508784
46 -0.40354429
47 -0.76409652
48 -0.28436510
49  0.08229647
50  0.63515448
51 -0.82162900
52  0.17124659

# los guardamos en formato tibble:
pca_scores <- as_tibble(pca_fit$ind$coord)
pca_scores # tabla con sus 52 observaciones (las 52 provincias) y las 6 comps como

# A tibble: 52 × 6
    Dim.1     Dim.2   Dim.3  Dim.4   Dim.5   Dim.6
    <dbl>     <dbl>   <dbl>  <dbl>   <dbl>   <dbl>
 1 -1.41   0.473     0.0962 -0.458  0.131  -0.379 
 2  3.38   0.540     1.92    2.42  -1.94    1.80  
 3 -0.617  2.61      0.208  -0.142 -0.0419  0.0504
 4 -1.44  -0.000844 -2.03   -0.113  0.468   0.0944
 5 -0.204 -1.95      1.30   -0.714 -0.369   0.0167
 6 -1.05   1.17      1.80   -0.854  0.0754 -0.575 
 7  1.53   0.260    -2.52    2.36   0.481  -0.567 
 8 13.7   -1.61     -0.867  -0.279  1.04    0.765 
 9  0.576 -1.51     -1.18   -0.383  0.798   0.466 
10 -1.20  -1.55     -0.892   0.989  0.602   0.532 
# … with 42 more rows

           # nuevas variables (columnas)

¿Cuál es la contribución de las variables originales en cada componente principal seleccionada?

pca_fit$var$contrib

                   Dim.1        Dim.2        Dim.3       Dim.4
Poblacion     8.61580869 6.216756e-04 2.494737e-01  0.28510690
Mortalidad    1.12979868 2.778556e+01 3.567076e+00  2.60231048
Natalidad     0.16505544 2.454404e+01 7.326830e+00  1.20275871
IPC           1.20805785 1.334604e+01 6.881053e+00 18.92662230
NumEmpresas   8.65432772 6.836923e-02 6.323813e-03  0.47800598
Industria     8.15804006 2.001441e-01 2.150098e-01  0.05499531
Construccion  8.60042932 2.059505e-01 1.459770e-02  0.06966501
CTH           8.57705398 1.127314e-02 2.385458e-01  0.07607058
Infor         7.92433313 1.765523e-01 4.178363e-01  4.90741583
AFS           8.55363151 2.707409e-02 1.564543e-01  0.84054061
APT           8.45209684 8.663137e-02 7.878255e-02  2.01245468
TasaActividad 1.30730170 1.093060e+01 1.315763e+01 21.44167757
TasaParo      0.01958975 2.130396e+01 1.500454e+01  4.84888573
Ocupados      8.66880711 2.835454e-02 5.901685e-04  0.36063054
PIB           8.46269279 1.307379e-01 1.405631e-01  1.78508842
CANE          0.03164342 9.251454e-01 4.313115e+01  7.75029908
TVF           8.50111166 2.558733e-04 9.952426e-01  0.19487506
VS            2.97022035 2.286851e-01 8.418293e+00 32.16259720
                     Dim.5        Dim.6
Poblacion     2.527731e-02  0.008577510
Mortalidad    5.830576e-02  0.113808680
Natalidad     1.459711e+01  9.825724095
IPC           3.099788e+01 12.301662219
NumEmpresas   4.053946e-03  0.004628499
Industria     2.251289e+00  0.379668213
Construccion  6.583274e-03  0.003921397
CTH           8.046371e-02  0.223154480
Infor         4.960879e-01  1.577980441
AFS           9.251174e-04  0.299366263
APT           1.006972e-01  0.640893842
TasaActividad 1.770063e+00 27.932330747
TasaParo      3.065031e-01  9.253743706
Ocupados      1.518964e-04  0.159122519
PIB           7.434893e-03  0.593372087
CANE          2.381808e+01 17.575176869
TVF           2.178721e-01  0.543501985
VS            2.526122e+01 18.563366447

¿Cuál de las variables es la que está peor explicada en función de las dos primeras componentes?

# Visualmente, viendo el gráfico del EJERCICIO 4 pto 3.a) (gráfico fviz_pca_var
# para Dim.1 y Dim.2) es inmediato responder que la peor explicada por comps 1 y 2
# es CANE.
 
# Numéricamente, la idea es también fijarnos en el cos2 total de ambas comps:
which.min(rowSums(pca_fit$var$cos2[,c(1,2)])) # el resultado es CANE

CANE 
  16 

Ejercicio 5:

Si tuviéramos que construir un índice que valorase de forma conjunta el desarrollo económico de una provincia, cómo se podría construir utilizando una combinación lineal de todas las variables. ¿Cuál sería el valor de dicho índice en Madrid? ¿Cual sería su valor en Melilla?

# Un índice es un único valor numérico que nos proporcionara una estimación del 
# desarrollo económico; es decir, que reflejara la mayor cantidad de información 
# posible contenida en las variables económicas, y que tomara un valor para cada
# distrito basado en los valores de esas variabs para cada provincia.

# Y hemos visto que la componente principal del análisis PCA hace precisamente eso,
# buscar una nueva variable (componente) que capture la mayor cantidad de información
# posible (en una única variable o un único valor, ie, un "índice").

# Por tanto, lo que debemos hacer para construir ese índice es aplicar el PCA a las 
# variabs orig de carácter económico.

# No obstante, el enunciado indica que usemos TODAS la variables, así que lo que
# haremos será utilizar el PCA que ya hemos realizado, y utilizar como índice
# económico la componente en la cual las variables económicas tienen más relevancia
# (ie, la componente que mejor capture la información económica, que es justamente 
# lo que queremos que refleje el índice).

# Para identificar esa componente, tendríamos que fijarnos en qué comp tiene los
# mayores valor total de cos2 (que recordemos nos daba una estimación de lo bien
# representada que está una variab orig en cada comp).
which.max( colSums( pca_fit$var$cos2[-c(1,2,3,17,18), ] ) ) # Dim.1

Dim.1 
    1 

# Además, en este caso en particular, y recordando el gráfico del EJERC 4 apartado 3.b,
# es inmediato ver que la comp que buscamos es la Dim.1, pues la mayoría de variabs
# de carácter económico tienen una elevada correlacion con la Dim.1.

# Por tanto, la Dim.1 constituye el índice económico que buscamos.Y el valor del índice
# para cada provincia será el valor en esa comp de cada provincia, esto es, el score:

aux <- which(provincias$Prov == "Madrid") # Madrid es la fila 30
pca_fit$ind$coord[aux, 1] # el valor del índice en Madrid es 16.78

[1] 16.7783

pca_fit$ind$coord[which(provincias$Prov == "Melilla"), ]

     Dim.1      Dim.2      Dim.3      Dim.4      Dim.5      Dim.6 
-2.2179459  4.7821957 -1.9051614 -2.2308254  0.7138309  1.4096604 

# el valor del índice en Melilla es -2.22

# [EXTRA] Podemos ver otras provincias:
pca_fit$ind$coord[which(provincias$Prov == "Barcelona"), 1] # 13.68

[1] 13.68269

pca_fit$ind$coord[which(provincias$Prov == "Valencia"), 1] # 4.77

[1] 4.769832

pca_fit$ind$coord[which(provincias$Prov == "Bizkaia"), 1] # 0.58

[1] 0.5764437

pca_fit$ind$coord[which(provincias$Prov == "Soria"), 1] # -2.40

[1] -2.399059

Ejercicio 6:

Calcula la matriz de distancias de los datos (realiza las modificaciones que sean necesarias antes). Representa dicha matriz con un mapa de calor. ¿Se detectan inicialmente grupos de provincias? ¿Cuántos?

# Empezamos el CLUSTERING ANALYSIS

# (Para visualizar mejor los gráficos en los apartados siguientes, vamos a 
# cambiar los id numéricos por las propias provincias)
rownames(provincias) <- provincias$Prov

# a) [EXTRA] mapa de calor de los datos (sin estandarizar) - no confundir con
# la matriz de distancias, que la representaremos después

heatmaply(provincias %>% select(-Prov),
          seriate = "mean", # Con seriate cambiamos el orden (OLO, mean...)
          row_dend_left = TRUE,
          plot_method = "plotly")

# Con los datos sin estandarizar, unos pocos datos de valor numéricos muy alto
# (concretamente del PIB, sobretodo Barcelona y Madrid) hacen que el resto de valores,
# en proporción, no se distingan visualmente. Por eso, luego lo veremos mejor en la
# representación con los datos estandarizados

# (este apartado continua después del gráfico)

# b) matriz de distancias (datos sin estandarizar)
d <- dist(provincias,
          method = "euclidean") # method indica la métrica

fviz_dist(d, show_labels = TRUE) + 
  ggtitle("Matriz de distancias (datos sin estandarizar")

# Ocurre que Barcelona y Madrid tiene valores muy diferentes al resto de provincias,
# y por eso sus distancias con el resto de ptos son muy grandes y, en esta escala,
# cuesta detectar diferencias entre las otras distancias. Los matices se verán mejor
# cuando representemos con datos estandarizados

# Aun así, sí se aprecian algunos grupos de ptos que están cerca entre sí y 
# distanciados del resto (clusters):
# - Madrid y Barcelona
# - Provincias (típicamente más "ricas"): Bizkaia, Sevilla, Navarra, etc.
# - Provincias (típ más "pobres"): Zamora, Teruel, Soria, Ourense, etc.
#  - Valencia es la única prov que no parece estar en ninguno de estos 3 grupos

# Para visualizar mejor los datos, los ESTANDARIZAMOS
provincias_std <- provincias %>% mutate(across(where(is.numeric), rescale))

# De nuevo, cambiamos los id numéricos por las propias provincias
rownames(provincias_std) <- provincias_std$Prov

# c) [EXTRA] mapa de calor de los datos (estandarizados)
heatmaply(provincias_std %>% select(-Prov),
          seriate = "mean", # Con seriate cambiamos el orden (OLO, mean...)
          row_dend_left = TRUE,
          plot_method = "plotly")

# con los datos estandarizar, ahora sí vemos muchos detalles que antes era 
# imposible detectar.Ahora ya podemos ver, por ej, por qué la primera distinción
# en grupo de indiv es por un lado Madrid y Barcelona, y por otro lado el resto.

# d) matriz de distancias (datos estandarizados)
d_std <- dist(provincias_std,
          method = "euclidean") # method indica la métrica

fviz_dist(d_std, show_labels = TRUE) +
  ggtitle("Matriz de distancias (datos estandarizados)")

# Observamos ahora una distribución de grupos más difusa que cuando vimos los datos 
# sin estandarizar, pero con la ventaja de que ahora podemos idenficar algunos grupos
# más; es decir, al estandarizar, hemos logrado mayor grado de distinción entre grupos:
# - Madrid y Barcelona siguen siendo un grupo alejado del resto
# - Valencia y Alicante son otro grupo
# - Ceuta y Melilla otro
# - Luego observamos dos grupos más, uno más grande (con subgrupos dentro) y el 
# otro más pequeño. Pero ahora vemos que las diferencias no están tan claras, incluso 
# parece que podría haber un grupo más (7 en total)... (lo mejor será usar los
# algoritmos de clustering)

Ejercicio 7:

Realiza al menos tres análisis de clúster jerárquico con distintos enlaces y comenta las diferencias. En cada caso visualiza el dendograma y comenta cuántos clústers recomendarías usar (haciendo uso de la información de apartados anteriores)

# algoritmos jerárquicos - función hclust

## 1) enlace SIMPLE o dist mín (vecinos más cercanos)
single_clust <- hclust(d_std, method = "single")

# - dendograma sin clusters
fviz_dend(single_clust) + 
  labs(title = "Dendograma (enlace SIMPLE)")

# - dendograma con k clusters

# Con la info que comentamos a raíz de la matriz de distancias representada, podríamos
# usar entre 3 y 6 clústers. (Para decidir mejor el número óptimo, podemos realizar
# luego un análisis para determinarlo, pero por ahora vamos a razonarlo observando
# el dendograma)

# Observando los cortes del dendograma y la info de apartados anteriores, pensamos
# que lo mejor sería usar k = 5. Veamos cómo quedaría:

col = c('darkgreen','green',"#E7B800","#2E9FDF",'purple','red','orange')

fviz_dend(single_clust, k = 5,
          cex = 0.8, # tamaño de texto del eje x
          k_colors = col,
          color_labels_by_k = TRUE, # color de texto del eje x
          ylim = c(-0.25,2),
          rect = TRUE) + # añade un rectángulo alrededor
  labs(title = "Dendograma (enlace SIMPLE)")

# Observamos que el método simple no ayuda a distinguir de forma clara grupos
# dentro del conjunto de provincias más allá de las que están claramente aisladas
# (Madrid, Barcelona, Valencia y Alicante). En el siguiente punto, veremos si otros
# métodos arrojan mejores resultados.

## 2) enlace COMPLETO o dist máx (indiv más lejanos)
complete_clust <- hclust(d_std, method = "complete")

# - dendograma sin clusters
fviz_dend(complete_clust) + 
  labs(title = "Dendograma (enlace COMPLETO)")

# - dendograma con k clusters
# Observando el dendograma, estimamos que, en este caso, k = 6 es lo óptimo

fviz_dend(complete_clust, k = 6, # probemos con 6 clusters
          cex = 0.8, # tamaño de texto del eje x
          k_colors = col,
          color_labels_by_k = TRUE, # color de texto del eje x
          ylim = c(-0.25,4),
          rect = TRUE) + # añade un rectángulo alrededor
  labs(title = "Dendograma (enlace COMPLETO)")

# Con el enlace completo sí observamos algo que nos recuerda más a lo que observamos
# con la matriz de distancias: por un lado, los 2 grupos "aislados" (Madrid y Barcelona,
# Valencia y Alicante), y por otro lado, los demás grupos, más cercanos entre sí,
# como si pudieran considerarse subgrupos de un mismo grupo. De hecho, podemos probar
# a reducir k (= 4 o 5) y no sería una mala división. Aunque, considerando que tenemos
# 52 observaciones con características bastante heterogéneas, nos parece más acertado 
# aumentar un poco el número de grupos (k = 6), para reflejar esa variedad.

# No subiríamos a 7 grupos ni más porque, como vemos en el dendograma, estaríamos 
# diferenciando algunos grupos pero dejando otros con distancias muy similares sin
# diferenciar, y eso supondría una división subóptima de los datos, en lo que a su
# interpretación se refiere (sería una mala categorización)

## 3) enlace medio o dist media entre observaciones (de distintos grupos)
average_clust <- hclust(d_std, method = "average")

# - dendograma sin clusters
fviz_dend(average_clust) + 
  labs(title = "Dendograma (enlace MEDIO)")

# - dendograma con k clusters
# Observando el dendograma, estimamos que, en este caso, k = 6 es lo óptimo

fviz_dend(average_clust, k = 8, # probemos con 6 clusters
          cex = 0.8, # tamaño de texto del eje x
          k_colors = col,
          color_labels_by_k = TRUE, # color de texto del eje x
          ylim = c(-0.25,4),
          rect = TRUE) + # añade un rectángulo alrededor
  labs(title = "Dendograma (enlace MEDIO)")

# en este caso, nos encontramos con un dilema, pues con k=8 justo dividimos los dos
# grupos de la derecha del gráfico, pero deja sin dividir, por muy poco, a Alicante
# y Valencia. Por un lado, ya mencionamos antes la importancia de no elegir un corte 
# que implique la división de algunos grupos pero no la de otros con distancias muy
# similares; no obstante, aquí el grupo de Valencia y Alicante está muy separado del 
# resto de observaciones, por lo que, aunque entre las propias Valencia y Alicante
# haya distancias similares a las de los otros grupos, creemos que tiene sentido
# (a nivel explicativo de los datos) cortar en k=8 y mantener esas dos provincias
# en un único grupo.

## Tras probar estos 3 métodos, lo que podemos concluir es que el método 2 es el que
## consigue distinguir los (6) grupos de forma más clara y diferenciada (el "punto"
## de corte es más elevado que en los otros dos métodos; su dendograma muestra una
## división más clara).

Ejercicio 8:

Echa un vistazo al final de las diapositivas de clase. ¿Qué número óptimo de clusters nos indican el criterio de Silhoutte y el de la varianza intracluster W? Detalla lo que consideres de ambos criterios y representa los individuos agrupados según el número de clusters elegido.

# Determinemos el número k óptimo de clusters (categorías) cuando no lo conocemos,
# empleando los criterios de silhouette y de varianza intracluster W ("wss")

# a) Silhouette

fviz_nbclust(provincias_std[ ,-1],
             kmeans,
             method = "wss") + # wss,
  geom_vline(xintercept = 3, # elegimos visualmente el k óptimo y lo pintamos
             linetype = 2) +
  geom_vline(xintercept = 5,
             linetype = 3) +
  theme_minimal() +
  labs(x = "nº clústeres (k)",
       y = "Variabilidad total intra-clústeres (W)",
       title = "Número óptimo de clusters basado en variabilidad total intra-clústeres")

# la línea vertical la hemos pintado a posteriori, tras determinar visualmente el
# k óptimo, escogiendo el k cuando la varianza intracluster W ya no se reduce de
# forma significativa ("Elbow method").

# En este caso en particular, hemos mostrado dos líneas porque consideramos que tanto
# k = 3 como k = 5 son candidatos a k óptimos. 

# b) Silhouette

fviz_nbclust(provincias_std[ ,-1],
             kmeans,
             method = "silhouette") + # wss,
  theme_minimal() +
  labs(x = "nº clústeres (k)",
       y = "índice Silhouette",
       title = "Número óptimo de clusters basado en el índice de Silhouette")

# la línea vertical la escoge y pinta la propia función. El método determina 
# k = 3 como nº óptimo. Este resultado coincide con el del método de enlace 
# SIMPLE visto anteriormente, y esto se debe a que el método silhouette se basa
# en la compacidad de los datos, y el método de enlace simple se basa en la distancia
# mínima entre vecinos, por lo que tiene sentido que arrojen conclusiones similares.

# Una vez detectados los dos posibles k óptimos (3 y 5), vamos a visualizar los
# datos agrupados tanto con k = 3 como con k = 5, para ver si podemos observar
# algo que nos ayude a decidir cuál de los valores es mejor, y también para ver cómo
# quedan distribuidos los clusters.

# Para representar, escogeremos el método de enlace COMPLETO, pues con lo visto en el
# ejercicio 7, consideramos que es el mejor método (de los 3 que hemos utilizado) de
# selección de clusters para nuestro dataset.

# k = 3
groups_3 <- cutree(complete_clust, k = 3)

fviz_cluster(list(data = provincias_std %>% select(-Prov),
                  cluster = groups_3), # groups
             palette = c('blue','red','green'),
             ellipse.type = "convex",
             repel = TRUE,
             show.clust.cent = FALSE,
             geom = "point") + # por dfto c("point","text"), pero text lo pondrá geom_text

  geom_text(label = rownames(provincias_std),
            aes(color = as.factor(groups_3)), #aes con as.factor para que se pinte bien
            check_overlap = TRUE,
            nudge_x = 0.25, nudge_y = 0.25,
            show.legend = FALSE) + # para que no pinte una "a" en la leyenda) +
  
  labs(title = "Cluster (complete)") + 
  theme_minimal()

# k = 5
groups_5 <- cutree(complete_clust, k = 5)

fviz_cluster(list(data = provincias_std %>% select(-Prov),
                  cluster = groups_5), # groups
             palette = c('blue','red','green','purple','orange'),
             ellipse.type = "convex",
             repel = TRUE,
             show.clust.cent = FALSE,
             geom = "point") + # por dfto c("point","text"), pero text lo pondrá geom_text

  geom_text(label = rownames(provincias_std),
            aes(color = as.factor(groups_5)), #aes con as.factor para que se pinte bien
            check_overlap = TRUE,
            nudge_x = 0.25, nudge_y = 0.25,
            show.legend = FALSE) + # para que no pinte una "a" en la leyenda) +
  
  labs(title = "Cluster (complete)") + 
  theme_minimal()

# Observamos lo que ya sabíamos, que el clustering con k = 3 solo separa 4 provincias
# (Madrid, Barcelona, Valencia y Alicante) del resto. Esto no nos parece una
# categorización demasiado útil, y además, vemos que k = 3 no captura para nada 
# la información de la comp Dim.2 (solo captura la de Dim.1), y esto es relevante
# porque recordemos que entre ambas, reflejan el 78% de la info de nuestros datos.

# Es por eso que consideramos más adecuado quedarnos con k = 5 (*)

# (*) si el método de clustering es un algoritmo jerárquiro (luego veremos lo que
# pasa con el algoritmo no jerárquico de k-means).

# Comentario: observamos que k = 5 da lugar a dos clusters que, vistos en la proyección
# de Dim.1 y Dim.2, se superponen en algunos puntos; cabe pues recordemos que 
# estamos viendo la proyección en 2 dims de puntos que están en un espacio n-dimensional
# (donde n es el nº de variabs orig), y que en ese espadio de mayor dimensión, los clusters 
# sí estarán separados.

Ejercicio 9:

Con el número de clusters decidido en el apartado anterior realizar un agrupamiento no jerárquico de k-medias. Representa los clusters formados en los planos de las Componentes principales. Interpreta los resultados y evalúa la calidad del análisis clúster. Explica las provincias que forman cada uno de los clusters y comentar cuales podrían ser las características socioeconómicas que las hacen pertenecer a dicho clúster

# ALGORITMO K-MEANS

# Probemos primero con k = 5:
kclust_5 <- kmeans(provincias_std %>% select(-Prov),
                 centers = 5, # k óptimo escogido
                 nstart = 25, # (*)
                 iter.max = 50) # iteracciones máximas del algoritmo

# representar los clusters en las Comp Princ:
fviz_cluster(list(data = provincias_std %>% select(-Prov),
                  cluster = kclust_5$cluster),
             palette = c("#2E9FDF","#E7B800","#00AFBB",'red','green'),
             ellipse.type = "convex", 
             repel = TRUE,
             show.clust.cent = FALSE,
             geom = "point") + # por dfto c("point","text"), pero text lo pondrá geom_text

  geom_text(label = rownames(provincias_std),
            aes(color = as.factor(kclust_5$cluster)), #aes con as.factor para que se pinte bien
            check_overlap = TRUE,
            nudge_x = 0.25, nudge_y = 0.25,
            show.legend = FALSE) + # para que no pinte una "a" en la leyenda
  
  labs(title = "Cluster (k-means) con k = 5") + 
  theme_minimal()

# (*) se ha incluido nstart = 25 para solucionar el siguiente problema:

# PROBLEMA: hemos detectado que, con k = 5, el análisis k-means nos devuelve 
# resultados muy diferentes cada vez que lo ejecutamos. En otras palabras, con k = 5,
# el algoritmo encuentra distintas formas de agrupar los datos en 5 clusters.

# Antes de descartar k = 5 como k óptimo por este motivo, recordemos que la función
# kmeans() utiliza centroides aleatorios (pues no se los estamos proporcionando).
# Podríamos "ajustar" un poco esos centroides para obtener resultados más
# consistentes.
# 
# Existen varias opciones, pero nos ceñiremos al temario de clase y
# nos limitaremos a cambiar el argumento nstart (por defecto nstart = 1) y lo
# fijaremos a nstart = 25, para que al incluir más centroides iniciales que los
# k clusters que queremos al final, veamos si conseguimos mejorar la consistencia.
# 
# Efectivamente, con nstart = 25 el resultado sí es consistente.

# Evaluemos la CALIDAD del análisis k-means (con k = 5)

# función silhouette:
sil <- silhouette(kclust_5$cluster, d_std) # kclust$cluster

row.names(sil) <- row.names(provincias_std) # filas, no cols! Del "1" al "150"

# Visualización
fviz_silhouette(sil, label = TRUE) +
  scale_fill_manual(values = c("#2E9FDF","#2edf6e","#E7B800","#FC4E07",'purple')) +
  scale_color_manual(values = c("#2E9FDF","#2edf6e","#E7B800","#FC4E07",'purple')) +
  theme_minimal() +
  labs(title = "Índice silhouette para k-means con k = 5") +
  theme(axis.text.x = element_text(angle = 90, # (girar etiquetas eje x)
                                   vjust = 0.5, hjust=1))

  cluster size ave.sil.width
1       1   17          0.22
2       2   13          0.36
3       3    2          0.61
4       4    4          0.26
5       5   16          0.21

# Explicación:
# La silueta mide la similitud entre una observación y las demás en su mismo cluster,
# en contraste con su similitud con las de otros clusters. Su valor oscila entre -1 y 1,
# de modo que un valor cercano a 1 indica que la observación está agrupada correctamente
# (ie, es muy similar al resto del cluster), mientras que cercana a -1 indica que 
# podríamos considerar cambiarla de cluster para mejorar el resultado.

# Tener muchos valores lo más cercanos a 1 posible nos indica que el
# clustering realizado es adecuado. Si por el contrario tenemos muchos valores bajos o negativos,
# podríamos considerar que el agrupamiento no es adecuado.

# En nuestro caso, para k = 5 vemos que todas las observaciones están por debajo de
# 0.5, excepto 2. Y muchas por debajo de 0.25. Luego veremos que este no es el mejor
# escenario, pues con otros k clusters se obtienen mejores resultados.

# En este punto, recordemos que escogimos k = 5 en lugar de k = 3 porque, para el
# método de clustering por enlace completo, solo separaba Madrid, Barcelona, Valencia
# y Alicante en 2 clusters, y todas las demás en un tercer y único cluster.
 # Pero vamos a probar el k-means con k = 3 y ver si el resultado difiere.

# k-means (k = 3)
kclust_3 <- kmeans(provincias_std %>% select(-Prov),
                 centers = 3,
                 # nstart = 25, # con k = 3 no ha hecho falta
                 iter.max = 50)

fviz_cluster(list(data = provincias_std %>% select(-Prov),
                  cluster = kclust_3$cluster),
             palette = c("#2E9FDF",'red','green'),
             ellipse.type = "convex", 
             repel = TRUE,
             show.clust.cent = FALSE,
             geom = "point") + 
  geom_text(label = rownames(provincias_std),
            aes(color = as.factor(kclust_3$cluster)),
            check_overlap = TRUE,
            nudge_x = 0.25, nudge_y = 0.25,
            show.legend = FALSE) + 
  labs(title = "Cluster (k-means) con k = 3") + 
  theme_minimal()

# CALIDAD (k = 3)
sil <- silhouette(kclust_3$cluster, d_std)
row.names(sil) <- row.names(provincias_std)

fviz_silhouette(sil, label = TRUE) +
  scale_fill_manual(values = c("#2E9FDF","#2edf6e","#E7B800","#FC4E07",'purple')) +
  scale_color_manual(values = c("#2E9FDF","#2edf6e","#E7B800","#FC4E07",'purple')) +
  theme_minimal() +
  labs(title = "Índice silhouette para k-means con k = 3") +
  theme(axis.text.x = element_text(angle = 90,vjust = 0.5,hjust=1))

  cluster size ave.sil.width
1       1    2          0.67
2       2   27          0.39
3       3   23          0.19

# Vemos que mejora ligeramente el resultado (fijarse en la línea roja punteada, que
# da la media).

# Sin embargo, sí hay algo que nos da que pensar:
# A diferencia de lo que ocurría con el método COMPLETO (algoritmo jerárquico), 
# cuando probamos con el algoritmo k-means (no jerárquico) con k = 3, el resultado
# sí refleja la diferenciación debida a la Dim.2. El problema es que ahora hay poca
# diferenciación debida a la componente principal, y eso no nos convence, ya que
# ésta es la comp que explica la mayor varianza de los datos.

# Por este motivo, vamos a probar también k-means con el k intermedio, k = 4.
kclust_4 <- kmeans(provincias_std %>% select(-Prov),
                 centers = 4,
                 nstart = 25,
                 iter.max = 50)
fviz_cluster(list(data = provincias_std %>% select(-Prov),
                  cluster = kclust_4$cluster),
             palette = c("#2E9FDF",'red','green','orange'),
             ellipse.type = "convex", 
             repel = TRUE,
             show.clust.cent = FALSE,
             geom = "point") +
  geom_text(label = rownames(provincias_std),
            aes(color = as.factor(kclust_4$cluster)),
            check_overlap = TRUE,
            nudge_x = 0.25, nudge_y = 0.25,
            show.legend = FALSE) + 
  labs(title = "Cluster (k-means) con k = 4") + 
  theme_minimal()

# Calidad (k = 4)
sil <- silhouette(kclust_4$cluster, d_std)
row.names(sil) <- row.names(provincias_std)

fviz_silhouette(sil, label = TRUE) +
  scale_fill_manual(values = c("#2E9FDF","#2edf6e","#E7B800","#FC4E07",'purple')) +
  scale_color_manual(values = c("#2E9FDF","#2edf6e","#E7B800","#FC4E07",'purple')) +
  theme_minimal() +
  labs(title = "Índice silhouette para k-means con k = 4") +
  theme(axis.text.x = element_text(angle = 90,vjust = 0.5,hjust=1))

  cluster size ave.sil.width
1       1    2          0.63
2       2   25          0.40
3       3    7          0.13
4       4   18          0.23

# El resultado es ligeramente mejor que con k = 5, y más o menos similar a k = 3.
# Por esto, tanto k = 3 como k = 4 serían k's óptimos desde el punto de vista de la calidad.

## NOTA: también se ha probado con k = 6 y k = 7, pero se obtenían resultados
# de calidad notablemnte peor, por lo que se han descartado como opciones.

# CONCLUSIÓN FINAL

# Con un análisis de clustering con el algoritmo de k-means, los mejores resultados 
# se obtienen para k = 3 y k = 4, siendo k = 3 el mejor, si tenemos en cuenta los análisis
# que vimos anteriormente basados en la varianza intracluster (wss) y en la silueta,
# como según el análisis de calidad (también basado en el criterio de la silueta) que
# hemos realizado en este apartado.

# No obstante, creemos que k = 4 también ofrece buenos resultados, y 
# podríamos considerar otros aspectos del problema o situación que estuviéramos
# estudiando, y decidir si nos conviene más tener 3 o 4 categorías.

Explica las provincias que forman cada uno de los clusters y comentar cuales podrían ser las características socioeconómicas que las hacen pertenecer a dicho clúster.

### Vamos a responder esta pregunta para el caso de k-means con k = 3.

# Previamente, veamos algunas cosas que podrían ser útiles en ciertas situaciones:

# Provincias del cluster 1:
provincias[which(provincias$cluster == 1),]

# A tibble: 0 × 19
# … with 19 variables: Prov <chr>, Poblacion <dbl>, Mortalidad <dbl>,
#   Natalidad <dbl>, IPC <dbl>, NumEmpresas <dbl>, Industria <dbl>,
#   Construccion <dbl>, CTH <dbl>, Infor <dbl>, AFS <dbl>, APT <dbl>,
#   TasaActividad <dbl>, TasaParo <dbl>, Ocupados <dbl>, PIB <dbl>,
#   CANE <dbl>, TVF <dbl>, VS <dbl>

# Provincias del cluster 2:
provincias[which(provincias$cluster == 2),]

# A tibble: 0 × 19
# … with 19 variables: Prov <chr>, Poblacion <dbl>, Mortalidad <dbl>,
#   Natalidad <dbl>, IPC <dbl>, NumEmpresas <dbl>, Industria <dbl>,
#   Construccion <dbl>, CTH <dbl>, Infor <dbl>, AFS <dbl>, APT <dbl>,
#   TasaActividad <dbl>, TasaParo <dbl>, Ocupados <dbl>, PIB <dbl>,
#   CANE <dbl>, TVF <dbl>, VS <dbl>

# Provincias del cluster 3:
provincias[which(provincias$cluster == 3),]

# A tibble: 0 × 19
# … with 19 variables: Prov <chr>, Poblacion <dbl>, Mortalidad <dbl>,
#   Natalidad <dbl>, IPC <dbl>, NumEmpresas <dbl>, Industria <dbl>,
#   Construccion <dbl>, CTH <dbl>, Infor <dbl>, AFS <dbl>, APT <dbl>,
#   TasaActividad <dbl>, TasaParo <dbl>, Ocupados <dbl>, PIB <dbl>,
#   CANE <dbl>, TVF <dbl>, VS <dbl>

kclust_3$centers # medias estandarizadas

  Poblacion Mortalidad Natalidad       IPC NumEmpresas  Industria
1 0.9269493  0.1926230 0.3182148 0.7531235  0.93315711 0.91207344
2 0.0673963  0.5697372 0.1702682 0.5021891  0.06293192 0.08728175
3 0.1295236  0.2566439 0.3120464 0.2956877  0.10264110 0.12689654
  Construccion        CTH      Infor        AFS        APT
1   0.92493092 0.99183626 0.79172581 0.85666694 0.85052570
2   0.08341758 0.07719223 0.02197145 0.05796761 0.04070624
3   0.10718001 0.13718569 0.03572563 0.10147636 0.06793674
  TasaActividad  TasaParo   Ocupados        PIB      CANE        TVF
1     0.7257989 0.1904479 0.90507944 0.86572133 0.1345357 0.94917422
2     0.4258563 0.1564716 0.06546121 0.05441053 0.2004695 0.09514499
3     0.5450924 0.6262214 0.10655128 0.07721645 0.3854134 0.16719253
         VS
1 0.4792530
2 0.1554777
3 0.2646706

kclust_3$cluster # clusters

 [1] 3 3 3 2 2 3 2 1 2 2 2 3 3 3 2 2 3 3 3 2 2 3 3 3 2 3 2 2 2 1 3 3 3
[34] 2 2 2 3 2 2 2 3 2 3 2 3 2 3 3 2 2 2 2

provincias$cluster <- kclust_3$cluster # esta col luego nos sirve para el aggregate


# MEDIA de cada variab orig para cada uno de los 3 clusters:
prov_agrupadas_media <- aggregate(provincias[ ,2:19], list(provincias$cluster), mean)

# RANGO de cada variab orig para cada uno de los 3 clusters:
prov_agrupadas_max <- aggregate(provincias[ ,2:19], list(provincias$cluster), max)
prov_agrupadas_min <- aggregate(provincias[ ,2:19], list(provincias$cluster), min)
prov_agrupadas_rango <- prov_agrupadas_max - prov_agrupadas_min # rangos
prov_agrupadas_rango$Group.1 <- c(1,2,3)


# ------------------------------------------------------------------------------
# PARTE 1: Explica las provincias que forman cada uno de los clusters

# Analizando los datos, hemos visto un fenómeno interesante: aparte del clúster de Madrid
# y Barcelona, los otros dos clústers separan las provincias entre mitad norte y mitad sur
# de España. Lo mostramos:
  
df <- provincias[,c(1,20)]
new_df <- df[order(df$cluster),]
new_df # esto nos muestra qué provincias pertenecen a cada cluster

# A tibble: 52 × 2
   Prov      cluster
   <chr>       <int>
 1 Barcelona       1
 2 Madrid          1
 3 Álava           2
 4 Asturias        2
 5 Balears         2
 6 Bizkaia         2
 7 Burgos          2
 8 Cantabria       2
 9 Coruña          2
10 Cuenca          2
# … with 42 more rows

# Vemos que efectivamente ocurre esa división entre NORTE y SUR, con 2 excepciones:
# - Guadalajara, que está en el cluster de las provincias del sur, pero geográficamente
#   está por encima de Cuenca (y Cuenca está en el clúster de las del norte - 
#   geográficamente está en el medio).
# - Tarragona, que está claramente en el norte, pero está en el cluster de las del sur.


# ------------------------------------------------------------------------------
# PARTE 2: comentar las características socioeconómicas que las hacen pertenecer a dicho clúster

# Con el siguiente gráfico, podemos comparar visualmente los clusters variable a variable:
# - mortalidad: la media del grupo 2 (norte) parece ser mayor que la dos grupo 3 (sur)
# - natalidad: la media del grupo 3 (sur) parece ser mayor que la dos grupo 2 (norte)
# Etc...
# La info que veríamos en estos gráficos aportaría, a fin de cuentas, la misma información
# que los valores de la media y el rango por clusters que vimos anteriormente, solo que
# aquí lo visualizaríamos y tendríamos más presentes los propios datos originales.

ggplot(provincias, aes(x = Prov, y = Natalidad, group = cluster)) +
  geom_point(aes(color = as.factor(cluster)), cex=3) + 
  labs(x = "Provincias") + 
  theme(legend.title=element_blank(),
        axis.text.x=element_text(angle=90, vjust=0.5, hjust=1))

# Así, podríamos analizar una a una cada variable (su gráfico, su media y su rango), 
# dictaminarqué variabs son (más) relevantes para clasificar las provincias en uno u
# otro cluster, en función de a qué otras provincias se asemejan.
# No obstante, lo más cómodo es usar PCA y representar en las 2 comps princ las provincias
# agrupadas en los 3 clusters. Lógicamente, esto es lo que veníamos haciendo anteriormente
# (precisamente porque es la manera óptima de representar, en 2Dim, toda la info posible).

# Esto es útil porque las comps están relacionadas con las variabs, como ya vimos en el
# EJERCICIO 4. De modo que, analizando el gráfico de PCA con los k = 3 clusters obtenidos
# con k-means, podremos identificar las características socioeconómicas que más contribuyen
# a diferenciar los clusters:

# k-means (k = 3)
fviz_cluster(list(data = provincias_std %>% select(-Prov),
                  cluster = kclust_3$cluster),
             palette = c("#2E9FDF",'red','green'),
             ellipse.type = "convex", 
             repel = TRUE,
             show.clust.cent = FALSE,
             geom = "point") + 
  geom_text(label = rownames(provincias_std),
            aes(color = as.factor(kclust_3$cluster)),
            check_overlap = TRUE,
            nudge_x = 0.25, nudge_y = 0.25,
            show.legend = FALSE) + 
  labs(title = "Cluster (k-means) con k = 3") + 
  theme_minimal()

# proyección de las variabs orig sobre las comp 1 y 2
fviz_pca_var(pca_fit, col.var = "cos2",
             gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"),
             repel = TRUE) +
  theme_minimal() + 
  labs(color= "Peso") +
  ggtitle("Coordenadas de las variables en las dos Comp Principales")

# -- CLUSTER 1:
# Este cluster está formado por provincias (Madrid y Barcelona) con:
# a) altos valores de Dim.1
# interpretar la comp 1 es más fácil que la comp 2: para una provincia (o cluster), 
# tener valores altos de Dim.1 se traduce tener valores altos de las variabs económicas
# que vimos que están altamente correlacionadas con la comp 1 (ie: Número de empresas, 
# Industria, Construcción, etc.). Lógicamente, Madrid y Barcelona tienen valores de
# estas variabs notablemente más elevados que todas las demás provincias.
# 
# b) valores positivos de Dim.2
# interpretar la comp 2 ya es más complicado, pues es combinación lineal de varias
# variabs (principalmente TasaParo, TasaActividad, Natalidad, Mortalidad e IPC) que
# no apuntan en la misma dirección ni sentido entre sí, y que no están tan altamente
# correlacionadas con ella como lo están por ej Industria y Construcción con la comp 1.
# 
# TasaParo, Natalidad y TasaActividad contribuyen a incrementar la Dim.2, mientras
# que Mortalidad e IPC contribuyen a disminuirla (IPC contribuye menos, ya que su
# correlación con Dim.2 no es tan elevada, y además su peso se reparte entre Dim.1 y 
# Dim.2).
# Veamos cuánto valen estas variabs en ambas provincias:
provincias[which( (provincias$Prov == "Madrid") | provincias$Prov == "Barcelona"),
           c("Prov","Mortalidad","Natalidad","IPC","TasaParo","TasaActividad")]

# A tibble: 2 × 6
  Prov      Mortalidad Natalidad   IPC TasaParo TasaActividad
  <chr>          <dbl>     <dbl> <dbl>    <dbl>         <dbl>
1 Barcelona       8.18      9.64  105.     17.2          61.8
2 Madrid          6.75     10.2   103.     16.3          63.9

# Veamos un ejemplo:
# Barcelona tiene Dim.2 mayor que Madrid, y sin embargo, también tiene una Mortalidad
# mayor. ¿Cómo es esto posible, si la mortalidad contribuye negativamente a la comp 2?
# (y lo mismo ocurre con TasaParo)
# Pues ocurre porque, recordemos, la Dim.2 no es solo la mortalidad, sino una combinación
# de variabs. Así pues, aunque Barcelona tenga mayor mortalidad, Madrid tiene menor IPC,
# mayor natalidad, y mayor TasaActividad.
# 
# En cualquier caso, al ser solo dos provincias en este cluster y además ambas estar
# relativamente cercanas en Dim.1 y Dim.2, podemos concluir lo siguiente: el cluster 1
# engloba a las provincias con valores elevados de las 11 variabs orig que están altamente
# correlacionadas con la comp 1, y que ya vimos que están ligadas al desarrollo económico.


# -- CLUSTERS 2 y 3:
# Por su parte, los clusters 2 y 3 tienen valores positivos de la comp 1, pero mucho menores
# que el cluster 1. Lo más relevante es que ambos clusters abarcan rangos similares en dicha
# comp 1 (despuntando un poco las provincias de Valencia y Alicante), lo cual refleja peores
# valores económicos que Madrid y Barcelona. Y también nos indica que para diferenciar o
# clasificar las provincias entre estos dos clusters, debemos centrarnos, ahora sí, en
# la Dim.2.
# 
# Las provincias del cluster 3 (sur) toman todas valores negativos de la comp 2
# (excepto Cáceres, pero igualmente está muy cercana a cero); mientras que las del cluster 2
# (norte) toman valores mayoritariamente positivos, aunque algunas están cercanas a cero
# (y algunas como Girona y Baleares son negativas).
# 
# A rasgos generales, ya hemos visto (con el ejemplo de Madrid y Barcelona) que el
# hecho de que ambos clusters se diferencien principalmente por la Dim.2 no nos ayuda
# a determinar cómo se traduce esa diferencia en términos de las variabs orig. Pero
# al menos sí podemos afirmar que esa diferencia se debe mayoritariamente a las variabs
# que más peso tienen en la Dim.2, a saber: TasaParo, TasaActividad, Natalidad, Mortalidad
# e IPC.

Ejercicio 10:

Comenta o concluye cualquier aspecto que consideres y que no se haya respondido anteriormente.

# Se han incluido/comentado algunas cosas adicionales ya en los ejercicios anteriores.
# De esta forma, se integraban en el orden seguido en la tarea.